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Abstract 

A detailed theoretical study of a heat-driven lamp has been performed. This lamp uses a 
plasma produced in a thermionic diode. The light is produced by the resonance transition of 
cesium. An important result of this study is that up to 30% of the input heat is predicted to be 
converted to light in this device. This is a major improvement over ordinary thermionic energy 
converters in which only ~1% is converted to resonance radiation. Efficiencies and optimum 
inter-electrode spacings have been found as a function of cathode temperature and the radiative 
escape factor. The theory developed explains the operating limits of the device. 
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I. Introduction 

The topic of this work is the direct conversion of heat to light. The method investigated is 
to use the cesium plasma of a thermionic discharge to produce resonance radiation. The 
motivation for this study is to determine the feasibility of the thermionic-photovoltaic energy 

converter proposed by D. L. Chubb. The reason cesium is used is that the frequency of its 

resonance radiation is closely matched to the band-gap of several semiconductors. This match 
permits efficient operation of photovoltaic cells in strong contrast to the inefficient operation of 
photovoltaics when powered by broadband solar radiation. 

The major result of this work is that up to 30% of the input heat to the discharge may be 

converted to resonance radiation even with very low escape factors for the radiation. A brief 

description of the technical w^rk is included in Appendix A. The details of the model, the 
calculations, and the results are given in Appendix B. 

II. Papers and Presentations 

• A paper on this work was presented at the 1983 IEEE International Conference on 
Plasma Science. The abstract of this paper is included as Appendix A. 

• F. Stefani completed his Master’s Thesis on this work. It is included as Appendix B. 

• A paper by F. Stefani and J. L. Lawless will be submitted to IEEE Journal of 
Plasma Science (in progress). 

III. Student Participation 

Student: F. Stefani 
Period: 1983-1983 
Degree: Master of Science 
Thesis: See Appendix B 
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The resonance radiation from a thermionic discharge can be 
end to generate electricity by employing photovoltaic cell* 
to collect this energy*. Since the Cb rmonnnee line la veil 
matched to the band-gap of many semiconductors, Bght-to- 
electridty eonv tml on efficiencies of 1096 or more are 
pomlble. Thb paper •will Investigate bow efficiently a 
thermionic discharge may he made to conve rt hast to light. 

Normally, an optimised thermionic e nerg y eoerve r ter 
produces very little naoannee radiation. Our theoretical 
studies show that at higher currents and longer inter- 
alectrode gaps, n significant portion of the input haul can be 
converted to line radiation even If the radiative eeeape factor 
la as low as 1%. Neglecting lend loams and optica) collection 
difficulties, efficiencies of ap to 1096 for the conve r sion of 
heat to radiation In the discharge are predicted. 

The t h eoretical studies modal the ele c tro n transport, km 
transport, and inelastic anllkioee leading to ionisation, 
raeemhinntiee, and emission of rwonance radiation. The 
transport at resonance radiation is treated eaiag an escape 
fleeter. The emission rate la eoopled to the loelaation and 
recombination rates through the Cb excited state kinetics. 
Diffusion and energy balances determine the electron 
temperature and deadly distribution. The mods! makes 
dear an unintuitive result that the maximum efficiency with 
which fine radiation la -—**•-** is only weakly depea dent on 
the optical escape factor. High sffidenelm are typically 
found at later- el ectrode distances on the order of 4 cm 

An experiment has been constructed to mea su re the 
efficiency for converting thermal energy to resonance 
radiation and abo the maximum discharge length. 
Measurements are In p ra gma. Available axperimeaul 
nauHe will be presented. 
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INTRODUCTION 


h has been suggested [ 1 ] that if line radiation from a thermionic device could be 
emitted with sufficient intensity, then a device employing photo cells to collect 
this energy could be made to generate electricity efficiently. This follows from 
the observation that the band gap energies of several common photo-voltaics 
are well matched to the line radiation energies of cesium diode converters. 


Resonance Radiation 

6 P 1 / 2 * 

1.39 eV 

6 P 3 / 2 - 68 ^ 

1.46 ev 

Excitation Band Gaps 

Silicon 

1.12 eV 

Indium Phosphide 

1.27 eV 

Gallium Aresenide 

1.35 eV 


Table 1 Radiation and band gap energies. 

This study investigates the efficiency with which resonance radiation can be 
generated from a thermionic device and identifies the operating conditions for 
which radiation output and the radiation efficiency are maximized. It has been 
tentatively established that radiation efficiencies on the order of 30% are 
possible even for optically dense plasmas. This result, if confirmed by 
experiment, means that the prospects for the proposed device are indeed 
promising. 

The first part of this report is a discussion of the model that was employed. This 
is followed by a brief outline of the solution procedure. The last part of this 
report is a presentation of the results and an attempt to relate them to the 
physical events occurring in the plasma. A fisting of the computer progran 
developed for this study is included in the appendix. 


AN OVERVIEW OF THE MODEL 

The work presented herein is based on an analytical model of a cesium plasma 
diode. This model, described In detail in reference 2. is summarized here 
primarily to review the assumptions made in the context of the problem and to 
provide some insights that help make the results more understandable. . 

A control volume is taken as the region of the Inter-electrode space extending 
from the top of the emitter sheath to the top of the collector sheath. Three 
governing equations are written for the plasma in the control volume. These 
are conservation equations for mass, momentum, and energy. A soluton 
procedure is employed that permits the optical thickness of the gas to be varied 
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as a parameter. Since the problem ie 1*D in the direction of current flow it is 
possible to reduce the governing equations to algebraic expressions. When 
these are solved simultaneously for the electron temperature, the electron 
density, and the operating current, a solution is obtained. 

From the onset, it was assumed that a thermionic device capable r! generating 
high levels of radiation would be operating in the ignited mode. This means 
that the ions and electrons in the plasma are created predominantly by electron- 
atom collisions in the control volume rather than by contact ionization at the 
walls. In this case emission effects are of relatively little importance. Of greater 
importance are the assumptions made about the plasma in the inter-electrode 
space. This work has assumed the following: 

1) The plasma is singly ionized and is composed of electrons, cesium 
atoms, and cesium ions. 

2) The plasma exists at two temperatures, one for the heavier particles and 
one for the electrons. 

3) Both temperatures are isothermally distributed in space. 

The two temperature plasma assumption is based cn the fact that most 
collisions are elastic, and very little energy is transferred from the fast moving 
electrons to the slower moving heavier particles. The isothermal assumption is 
in effect a statement that the thermal conductivity of the plasma is high. This 
assumption has been confirmed experimentally as well as through numerical 
tests [3]. 


AMBIPOLAR DIFFUSION EQUATION 

The Ambipolar Diffusion Eqn. is a mass balance for the electrons and the ions 
in the control volume. It is derived by writing the conditions for mass 
conservation for the steady state plasma in terms of Ionization and 
recombination rate constants S, and a. 

v*r, * v»r, * S Ng t n - an 3 (1 ) 

To reflect conditions in the device transport equations for the electrons and ions 
must include the effects of electric fields as well as diffusion effects. In one 
dimension they are written as follows: 



r# n cTF/dx + (kT./e) dn/dx) 
r i ■ii|(-n (fF/dx ♦ (kT./e) dn/dx) 


( 3 ) 
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where: 


r.,r, 

e 

n 


¥ 

k 

T.T, 


electron and ion fluxes 
electron and ion mobilities 
elementary charge 
electron density - ion density 
electron motive 
Boltzmann's constant 
electron and ion temperatures 


Substituting the transport equation for the ions (eqn. 3) into the conservation 
equation (eqn. 1) yields a diffusion equation for the ions. 


- d/dx((kT,/e) dn/dx - n d47dx) * S N g , n ■ an 3 (4) 

Substuting the electron transport (eqn. 2) into this expression eliminates the 
electric field term and results in a conservation equation in terms of the coupled 
motion of electrons and ions: 


-D a d 2 n/dx 2 +r,d(p/p,)/dx-«= (1 + p/p.,)[S N fli n - an 3 ] (5) 

where D. - ^ ( kT, + kT, ) / e (6) 

This equation is called the Ambipolar diffusion equation and D a is the ambipolar 
diffusivity. The ambipolar diffusion equation can be simplified by neglecting the 
ratio of the ion and electron mobilities ji/p, (since it is much less than one) and 
by making the reasonable assumption that this ratio does not vary much across 
the inter-electrode distance. Assuming also that D a is constant reduces the 
ambipolar diffusion equation to the following form: 

-D a d 2 n/dx 2 * S N fll n • an 3 (7) 

Solving the ambipolar diffusion equation subject to the appropriate boundary 
conditions (Appendix B) yields the following solution and associated 
eigenvalue condition: 

n(x) ■ n (max) sn(2KX/d) (8) 

SN 0t - (4K 2 )D a /d 2 ♦ (1/2) a n 2 (9) 

where n (max) is the maximum electron density. 

Beth these equations are of interest and reveal different things about the 
plasma. The solution to the ambipolar diffusion equation gives the distribution 
of the electrons in the inter-electrode space. The eigenvalue condition gives 
the correspondence between electron temperature and the maximum density 


for which ionization and recombination are euch to balance the loss of particles 
to e walls by diffusion. This correspondence is termed the ignited condition. 

The function sn that appears in the solution of the ambipolar diffusion equation 
is the Jacobi elliptic function [4]. It is characterized by the quarterperiod K, 

whose value can range from x/2 to infinity. Both limits represent important 
cases in the problem we are solving. In the limit where K- x/2, the sn function is 
exactly the sin function and the electron distribution is similar to that shown in 
fig. 1. This case corresponds tu conditions where recombination is negligible, 
usually at lower values of the electron density. The other limiting case is that 
where the electron density is very close to the Saha value. This condition is 
referred to as the Saha limit. In this case the electron distribution departs 
significantly from the sin profile, being much steeper at the boundaries and flat 
at the center, as in fig. 2. 

In non-Saha conditions the quarter period is calculated from the ratio of the 
electron density to the Saha density; first by calculating the modulus m, and 
then by taking the arithmetic geometric mean to find K [5]. 

m ■ 1 / ( 2 S N fli /( a n 2 ) - 1 ) (10) 

In the Saha limit, the modulus approaches an asymptotic value of unity and 
round-off error in the a.g.m. calculation of K becomes appreciable. In this case 
K is calculated from the following expression: 

K-nd(<x/8mD a ) 1 ' 2 , m=>1 (11) 

This is derived by solving eqn. 10 for SNg,, substituting this expression into the 
eigenvalue condition (eqn. 9) and solving for K. 


TRANSPORT EQUATION 

The second governing equation is a momentum balance for the electrons. It can 
be written as a generalized Ohm's law and is referred to in this work as the 
transport equation. The transport equation is obtained by solving the electron 
flux (eqn. 2) for the electric field and then integrating this over the length of the 
inter-electrode distance to obtain the voltage drop across the plasma. 

r, - - n dY/dx + (kT^e) dn/dx) (12) 

d'F/dx - - iVp. - (kT,/e) (1/n) dn/dx (13) 

T 0 -Ti--Jr, (1/p.n)dx - (kVe) In (n^) 


( 14 ) 


(15) 


¥<> • *i - JR- (kT«/e) In (n,/n 0 ) 

where n 0 and a* art tha alactron densities at x ■ 0 and x ■ d. 

The first term in tha pracading expression is tha voltage drop across tha plasma 
due to ohmic heating. Tha second term represents the voltage drop due to 
concentration gradients (diffusion drop). To calculate the resistance in the 
plasma R, the first term in eqn. 14 must be integrated over the length of inter- 
electrode space. If the resistance is to be valid over a broad range of conditions 
then the electron mobility that appears in the integral must include the effects of 
both electron-atom collisions and electron-ion collisions. 

In evaluating this integral it is convenient to make use of an approximate 
superposition rule for mobilities: 

1/p. -1/p„ +1/11.1 (16) 

where p.. and are the mobilities resulting from electron-atom 
collisions and from electron-ion collisions. 

Thus: R-R^ + R .1 

R ■ J 1/(ep.n)dx - J 1/(e|i,.n)dx + J 1/(ep.in)dx (17) 


Since conductivity due to electron-ion effects (o #i - eji^n) is a function of T. and 
independent of the electron density it may be approximated as a constant and 
removed from the integration. For the conditions of interest to us the resistance 
from electron ion effects may be calculated from the following expression [6]. 

R*i - 

where o.i - .15085 T.^/InjA) (18) 

and A- 12389 V*(n/2)* 1/2 

The electron-atom mobility is also a function of temperature and independent of 
the electron density so it too may be taken out of the integration. It is: 

p M - vXe/3KT (19) 

In the above expression the mean velocity of the atoms is token from kinetic 
theory as v-V(8KT/s n) and the mean free path is based on cross-sections 
found in the literature. This gives the following expression for the resistance: 


( 20 ) 


R « i/ji # J(1/ n ) d x + d/o #i 

If the appropriate substitutions are made the last remaining integral can be 
evaluated and the following expression is obtained. 

R - d/(2n, a n) In [4n* / n 0 n d (1 -m)] ♦ d/o ti (21 ) 

*^he preceding discussion shows that the electron-atom resistance is inversely 
proportional to the electron density while the electron-ion resistance is 
practically independent of it. This being the case, at high electron densities 
resistance depends only on temperature. The general behavior of resistance 
versus electron density is plotted in figure 3. This behavior accounts for some 
important trends that appear in the results. 

Often times the voltage across the electrodes is a more meaningful parameter 
than the voltage drop across the plasma. For this reason it ic more useful to 
express the transport equation in terms of the former quantity. This requires that 
the voltage drops across the electrode sheaths be added to the voltage drop 
across the plasma. It is assumed .that positive (electron retaining) sheaths are 
formed at the electrodes and that the flux at the sheaths is equal to the random 
thermal flux times a Boltzmann factor. 

coll: J - enjv^ e (** v * /kT *) => V, - kVe ln(en d v^4J) (22) 

emit: J « J E - en d v./4 e (■•v«/kT«) => v k «kT./e In [en d v^4(J E -J)] (23) 

When the sheath voltages are added to the expression for the plasma drop 
derived above several of the terms combine ahd the result is the transport 
equation in its final form: 

V d « JR + (kVe) In [J/ (J E -J)] (24) 

The first term still represents the ohmic drop while the second term now 
incorporates the combir.od effects of the sheath and diffusion drops. 


ENERGY EQUATION 

The third necessary equation is an energy balance for the system. In its final 
form this equation states that the electrical power loss in the plasma is caused 
by three separate processes. Energy is lost as thermionically emitted electrons 
deposit onto the collector, as electrons and ions that are created in the plasma 
diffuse to the walls and lastly as radiation out: 


V d J * 2KJ E (T # -T E )/e ♦ P dW + Prad 


(25) 
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The fate of thermionically emitted electrons is illustrated in fig 4. These enter 
the control volume at the emitter temperature and at a rate corresponding to the 
saturation current given by the Richardson-Dushmann equation. Neglecting 
back-emission from the collector, the operating current is equal to the emitted 
current minus the backscattered current. Each of these currents has an energy 
flux associated with it consisting of potential and kinetic energy terms. The net 
sum of these fluxes is the first term in the energy equation. This term is only a 
function of the electron temperature and is referred to as the kinetic energy term. 

The sum of diffusion and radiation losses is the power required to maintain the 
plasma ignited. This is calculated from the net rate of ionization established by 
the the ambipolar diffusion equation and can oe expressed as the sum of the 
terms and P rad . 


p diff « 4 E 0 K D, N,/d where E 0 is the ionization energy (26) 

Pr.d-EepNepgA (27) 

As one might expect the diffusion losses are inversely proportional to the inter- 
electrode distance. The power lost by radiation is equal to the product of the 
number of atoms in the first excited level, the frequency at which they emit and 
the energy of this radiation: 

The escape factor g, is used to modify the Einstein A-coefficienf so that (for 
example) if the gas is optically dense, g is close to zero and very little energy is 
lost in the form of radiation. To calculate the total number of atoms in the 6p 
state the expression for the local 6p population (eqn. 12) is integrated over the 
control volume. This integral cannot be evaluated analytically and so it is 
necessary to perform the integration numerically. In contrast to diffusion losses 
radiation losses are proportional to the inter-electrode distance. Also it should 
be understood that P rad represents radiation emitted in all directions not just the 
radiation which is available for collection by photo-panels. 


SOLUTION PROCEDURE 

In summary of the proceeding section.', the model reduces to three nonlinear, 
algebraic equations in three unknowns. These are: 

The ambipolar diffusion equation. 

D a d 2 n/dx 2 = SN g ,n - an 3 (7) 

The transport equation. 

V d - JR + (kT,/e) in [J/ (J E -J)] 

And the energy equation. 


( 24 ) 


1 

< 
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V d - (J/J E )2K(T.-T E )/e + P diff /J E + P^Je (25) 

A solution of model requires that these be solved simultaneously. Although any 
method of solution will do, the method outlined below is the one that yields the 
most insight into the problem. This method takes advantage of the fact that the 
ambipolar diffusion equation is a function of only two of the three variables (T e 
and n) and as such can be solved first for the ignited condition which the other 
two equations must satisfy (fig. 5). This reduces the problem to that of solving 
two equations in two unknowns. These can be solved simultaneously, however 
It is useful to solve them independently for the variable J and then to search for 
solution points (fig 6). This allows one to better understand why solutions exist 
where they do. 


RESULTS 

In the study the escape factor, the emitter temperature, and the inter-electrode 
distance were varied as parameters while other conditions were fixed at the 
following values. 

Ion Temperature * 1500 K 
Emitter work function « 2.4 V 
Plasma drop > 1 V 
Electrode area (planar) « 1 cm 2 
Number density of plasma * 10 16 /cm 3 

Briefly, the rationale that went into this selection was the following. 1 ) The 
escape factor is necessarily a parameter since the optical density is 
independent of the macroscopic quantities calculated in the model. For this 
reason it is also true that some of the curves, especially those for higher values 
of the escape factor, may not correspond to physically achievable solutions. 2) 
Varying the emitter temperature is equivalent to varying to power of the device. 
At the onset it was not clear what the relationship between power and efficiency 
would be. 3) Finally, the inter-electrode distance has a direct effect on the area 
available for radiation. Although distance does not appear explicitly in the 
plotted results it is perhaps the most important variable in determining the 
radiation output. One important result is that seme optimal solutions involve 
inter-electrode spacings which are significantly greater than normally 
encountered. 

Regarding the other conditions, the only one that warrants any mention is the 
plasma drop of one volt. It is this characteristic that sets the proposed device 
apart from conventional thermionic generators and is largely responsible for the 
higher radiaton eficiencies . This is discussed further in the next section. 

The results of the study are summarized in figures 7,8, and 9. Three degrees of 
optical density were investigated corresponding to escape factor values of g«1 , 
g-.l , g«.0l . For all three cases the emitter temperature was varied from 1500K 
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to 2000K. And for each combination of g and T E the inter-electrode distance 
was varied from the minimum value for which a solution existed to greatest 
possible value for which a solution existed. Since there is no simple criteria for 
determining beforehand where solutions exist the method used was in effect 
trial and error, letting the computer determine where solutions were 
impermissible. An important part of understanding the results involves 
understanding why under certain conditions solutions fail to exist. 

in each case studied the three governing equations were solved for the electron 
temperature, electron density and operating current by the method outlined 
above. From this information the radiation output and the efficiencies were 
calculated. Radiation output is defined as P rad .and the overall efficiency is 
defined as the ratio of the radiation output to the total power consumed by the 
device. 

*1 ■ Prw/l p r*j + p difi + J ( ♦a + 2kT # /e )] (27) 

where $ A is * he work functon of the anode. 

The volumetric efficiency was also calculated. This is a figure of merit defined 
as the ratio of the radiation output to the power dissipated inside the plasma. 


T)vc - P«*/JV d (28) 

The data points that are plotted correspond to the solutions for which the 
efficiency was the greatest, in each case these coincide with the solutions for 
which the interelectrode spacing was also the greatest. That this is hardly a 
coincidence is taken up in the following section. 


DISCUSSION 

Two results of this study are of particular interest. These are that radiation 
efficiencies on the order of thirty percent are attainable and that this maximum is 
independent of the optical density of the gas. High radiation efficiencies result 
from the greater plasma drop that characterizes the proposed device. A greater 
plasma drop means that more energy is dissipated inside the plasma, 
additionally it happens that a larger fraction of this energy is dissipated as 
radiation. That this should be the case is not obvious and is explained in the 
following section. 

One important consequence of a higher plasma drop is that solutions exist at 
greater inter-electrode spadngs. This is because larger plasma dro>x can 
overcome a greater resistance. Ultimately it is the larger inter-e tetrode 
spacings that account for the higher radiation efficiencies. As mentioned 
earlier energy is lost from the plasma in three seperate processes. These are 
illustrated in figure 10 which shows the contributions of the various terms to the 
energy equation. The kinetic energy term describes the energy that is lost as 
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emitted electrons are deposited onto the collector. P djff represents diffusion 
losses to the walls and P rad is the energy lost as radiation. Since the electron 
temperature is determined by the ignited condition, the kinetic-energy term is 
independent of the plasma drop. Radiation and the diffusion losses on the 
other hand depend strongly on the inter-electrode spacing. At short distances 
diffusion losses are very high. At large distances radiation losses dominate. 
From this discussion one can understand why at higher voltage drops a greater 
percentage of the plasma losses will occur as radiation and consequently 
radiation efficiencies will be greater. 

This does not yet explain why the maximum radiation efficiency is independent 
of the optical thickness of the gas, for Intuitively, one would expect optically 
dense gasses to be less efficient at producing radiation. The explanation lies 
two parts. First, optically dense solutions exist at greater distances then 
optically thin solutions. This can be understood by referring once again to figure 
10 and noting that since Prad is proportional to the escape factor g, a smaller 
slope on the Prad line will cause the energy equation to intersect the tranport 
equation at a greater distance. The consequence of this is that net radiation 
output is not proportional to the escape factor. In some cases studied, reducing 
the escape factor by a factor of 100 had the effect of reducing the radiation 
output by only a factor of 5. In addition to this there is the fact that optically 
dense gasses produce their maximum radiation at lower value of the emitted 
current (fig 9) and hence at lower input powers. 

It is interesting to compare the proposed device with a Carnot engine operating 
between the two electrode temperatures. If we take the collector temperature to 
be such that the emitted current from the collector is one percent of the emitted 
current from the emitter, then the device is operating between 1500K and 745K. 
The Camot efficiency between these two temperatures is 50.3% and so our 
device is operating at 56% of Camot efficiency. 

The results then indicate that high radiation efficiencies can be obtained from 
cesium thermionic devices provided that there is an adequate votage drop 
across the plasma. The value of one volt that was used in the study is a 
reasonable value in that it can be generated without any electrical power input 
to the device. High efficiencies come by virtue of the higher inter-electrode 
distances which the solutions will accomodate. One particularly positive result 
is that efficiencies on the order of 30% can be obtained even with optically 
dense gasses. Since this is most likely to be the case under conditions of 
interest, this study lends some theoretical support to the proposed device. 
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EFFICIENCY VS. EMITTER TEMPERATURE 



1500 1600 1700 1800 1900 2000 

EMITTER TEMPERATURE K ' 

Figure 7. Efficiency vs. emitter temperature. 


VOLUMETRIC EFFICIENCY VS. EMITTER TEMPERATURE 
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Figure 8. Volumetric efficiency vs. emitter temperature. 
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Figure 1 1 . Maximum inter-electrode spacing vs. emitter temperature. 


APPENDIX A: REACTION RATE CONSTANTS FOR EXCITED STATE 

KINETICS 


To derive expressions for the ionizaton and recombination rate coefficients S 
and a, it is necessary to know the atomic states that are present in the plasma 
and the reactions that occur between these states. In this work three states and 
three reactions (shown below) were considered important in determining the 
ionization/recombination kinetics of the piasna. 


Cs 6 , 4 e* 

«=> 

CS 6 p 

4 

e* 

Csgp 

S» 

Cs 6 , 

4 

hv 

Csep ♦ e- 

« 

Cs* 

4 

2e- 


Based on these reactions the rate equations for the atomic states are written as 
follows: 


dN 6# /dt ■ 'Ng t nK£ t .gp 4 N^nK^, 4 N^gA 


dNgp/dt ■ -cN 6 p K € p^, ♦ nN(|K(| 4 p * N^gA - nNgpKgp_ c 4 n^Ko^p 

where K 6 ,. 6p is the foreward rate constant for the reaction between 
the ground state and the 6p state, and Kt P .e B the backwards rate 
constant. The subscript c refers to the continuum. 

Under quasi-steady conditions the populations of the excited states change 
slowly in time. With this assumption the rate equation for the 6P state reduces 
to an algebraic expression that can be solved for the 6P population: 

N6p * Ro ♦RiNg, 

where : R© * / (nK^^ 4 riK^ 4 g A) 


R t -nN^Ke^sp / (nKgp^t 4 4 gA] 


Substituting the expression for N* p Into the rate equation above then yields the 
desired form from which the coefficients S and a.can be factored: 


dN 6# /dt * — nSNgg 4 an 3 


where: 


S ■ K||^p -R( ( ♦ flA / n ) 

o - R® ( nK^ ♦ gA ) / n* 

To evaluate these expressions it is necessary to know the troreward and 
backward rate constants for the reactions between adjacent states. 

K 6$ . 6p is the rate constant for collision*! excitation from the ground state to first 
excited level. It is found by integrating the product of the cross-section for the 
reaction [ref. 7] and the velocity distribution of the electrons (assumed 
Maxwellian) over the range cf energies. Assuming that the cross-section 
increases linearly from the threshold the result is: 

Ke. ep - exp (-C V 0 2 ) [B (V 0 « ♦ V 0 */C ♦ 1/C 2 ) + A (V 0 2 ♦ 1/C )]/2C 

where A « (m^ftkT.J^TtOo 

B - (m 0 /2JCkT # ) M 2Jtm # Aa/AE 

C - m*/2kT # 

and V 0 -V2m # /E, E-1.33 eV 

O 0 « -1.02x1 0‘ 14 cm 2 


Ao/AE * 71X10* 16 cm 2 /eV 

From this expression it is possible to obtain an expression for the backward 
reaction, K 6p4t . The quasi-equilibrium assumption states that: 


If equilibrium with the continuum is also assumed then the ratio of the 
populations can be gotten from the Saha equation: 

NesMsp .(1/3) exp(Ee.VkT.) 

This gives the following expression for K*p_e,: 

K«p-«* " *W«p 0*3) WX P(E«»-*AT«)- 

Kc^p is an effective rate constant for collisions! de-excitation from tho continuum 
to the 6P state. This is given in reference 8 as: 


Kc-6p ■ - 4t > (e 2 /(4jt€o^T,)) 5 (Skl./jim,) 1 ® 


Invoking once again quasi-equilibrium and the Saha equation yields: 
»Vc - Kc* P (2nh 2 /mkT,) 3/2 exp (Esp/kT.^Up + 1 ) 


APPENDIX B: BOUNDARY CONDITIONS 


In this work two different sets of boundary conditions are used for the electron 
den**<ty at the electrodes. In the lowest order approximation, the walls are 
considered absorbing so that: 


n 0 - 0 n d « 0 


where the subscripts 0 and d are used to denote the emitter are collector 
sheaths. 

This assumption works well in finding the solution to the ambipolar diffusion 
equation. It is not adequate however for expressions in which terms on the 
order of X/d are important. This is the case when evaluating the resistance of 
the plasma. For here the assumption of zero electron density at the walls will 
lead to an infinite resistance. 

Appropriate nonzero boundary conditions are derived in reference 2 by 
developing expressions for electron and ion fluxes at the walls. These are the 
following: 

^ - 4fWj Vr/T. (D a dn/dx + (p/pjr # ) 


n d - 4fWj VT/T. ( -D.dn/dx - (n/p.)r.) 


Where: fl is the unit vector normal to the electrodes. On the camode surface fi is 
directed into the plasma. On the anode it is directed out of the plasma. T, is the 
electron flux emitted thermionically from the cathode, given by the Richardson- 
Dushmann equation 
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PROGRAM PL1 

WRITTEN BY F. STEFANI, 1984. 

THIS PROGRAM CALCULATES THE OPERATING CONDITIONS FOR A THERMIONIC 
DEVICE BASED ON AN ANALYTICAL MODEL DEVELOPED BY J.L. LAWLESS. 
SOLUTIONS ARE FOUND BY FIRST CALCULATING THE TEMPERATURE THAT 
SATISFIES THE IGNITED CONDITION AT A GIVEN ELECTRON DENSITY, AND 
THEN SOLVING THE TRANSPORT AND ENERGY EQUATIONS FOR THE OPERATING 
CURRENT BASED ON THAT TEMPERATURE. THIS IS DONE FOR A RANGE OF 
ELECTRON DENSITIES AND THE POINT WHERE THE THREE EQUATIONS ARE 
SATISFIED SIMULTANEOUSLY CAN BE FOUND BY INSPECTION OF THE OUTPUT. 
THE PROGRAM ALSO CALCULATES THE RADIATION EFICIENCY OF THE DEVICE. 

[REF: AN ANALYTICAL MODEL OF THERMIONIC DISCHARGES, J.L. LAWLESS, 

CARNEGI E-MELLON UNIVERSITY] 

SOME ADDITIONAL COMMENTS: 

• THE FOLLOWING QUANTITIES ARE ASSIGNED FIXED VALUES IN THE PROGRAM: 


PLASMA VOLTAGE DROP VD 
ION TEMPERATURE TI 
# OF NUCLEI HNUC 
ELECTRODE AREA AELEC 


1.0 [V] 

1500. [K] 

1.0E16 [CM*— 3] 

1.0 [ CM * * 2 ] 

(WFC; CATHODE WK. FTN. ) 


* ANODE WORK FUNCTION: WFA = WFC - VD 


• RECOMBINATION IS ASSUMED THROUGHOUT. 

• N6P IS COMPUTED BY NUMERICALLY INTEGRATING THE LOCAL 6P POPULATION 
OVER THE INTER-ELECTRODE SPACE. 


RADIATION EFFICIENCY IS 


CALCULATED IN TWO WAYS: 

EFF1=PRAD/ (PRAD+PDIFRJ (VA+2KTE/CHG) ) 
EFF2 =FRAD/ ( RJ * VD ) 


• SIMPSONS RULE HAS BEEN MODIFIED TO TAKE ADVANTAGE OF SYMMETRY. 

* OUTPUT IS WRITTEN TO FILE PL1.DAT 


C 

C 

n 

v* 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C LIST OF VARIABLES 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


INPUTS: 

DIST 

TEMIT 

WFC 

G 

OUTPUTS : 

NE 

XN1 

TE 

SMI 

XJEN 


INTER-ELECTRODE GAP 
EMITTER TEMPERATURE 
EMITTER WORK FUNCTION 
RADIATION ESCAPE FACTOR 


ELECTRON DENSITY 
ELECTRON DENSITY AT THE ANODE 
ELECTRON TEMPERATURE 
1 - MODULUS 

OPERATING CURRENT BASED ON ENERGY EQN. 


[CM] 

[K] 

[V] 


[CM* — 3 ] 

[ CM * * - 3 ] 
[K] 

[A/CM*»2] 


c 

RJ 

OPERATING CURENT BASED ON TRANSPORT EQN . 

[A/CM*«2] 

c 

R 

PLASMA RESISTANCE 

[OHMS] 

c 

REA 

RESISTANCE DUE TO ELECTRON-ATOM COLLISIONS 

[OHMS] 

c 

RE I 

RESISTANCE DUE TO ELECTRON-ION COLLISIONS 

[OHMS] 

c 

VK 

CATHODE SHEATH VOLTAGE 

[V] 

c 

VA 

ANODE SHEATH VOLTAGE 

[V] 

c 

EFF1 

RADIATION EFFICIENCY (SEE HEADER) 


c 

EFF2 

VOLUMETRIC EFFICIENCY " 


c 

SAHA 

FLAG FOR SAHA CONDITION: 


c 


SAHA=0 => PLASMA IS RECOMBINING 


c 


SAHA=1 => RECOMBINATION NEGLIGIBLE 


c 


SAHA=2 => " " " 

SMI = 0 

c 

c 

PHYSICAL 

CONSTANTS: 


c 

c 

ME 

ELECTRON MASS 

[KG/1.602E-15] 

c 

MI 

ION MASS 

[KG/1.602E-15] 

c 

CHG 

FUNDAMENTAL CHARGE 

[C/1.602E-19] 

c 

K 

BOTZMANNS CONST 

[eV/K] 

c 

ESP 

ENERGY Cs6p-Cs6s 

[ev] 

c 

ESP 

ENERGY Cs6s-Cs6p 

[eV] 

c 

ACOEFF 

EINSTEIN A COEFFICIENT 

[1/S] 

c 

LS 

DEGENERACY OF THE GROUND STATE 


c 

LP 

DEGENERACY OF THE FIRST EXCITED LEVEL 


c 

SIGIA 

CROSS-SECTION FOR ION-ATOM COLLISIONS 

[CM**2] 

c 

SI GEA 

CROSS-SECTION FOR ELECTRON-ATOM COLLISIONS 

[CM**2] 

c 

SIGO 

Y- INTERCEPT OF CROSS-SECTION VS. ENERGY CRV 

. [CM«»2] 

c 

DLTXCT 

SLOPE OF CROSS-SECTION VS. ENERGY CURVE. 

[CM**2/EV] 


c 

EXTERNAL F 
EXTERNAL SIMPER 
EXTERNAL THIN 
EXTERNAL XJT 
EXTERNAL TDIFF 

REAL KSP , KPS , KSPA , KSPB , KCP , KPC , NGS , NE , N6P , MUI , MUE , MFPI , MFPE 
REAL M I i KSPCON , ME , HD AP , PI , LS , LP , K , NNUC , JE , JNET , QP , MOD ( JION 
DIMENSION DDIST(15) ,DTEMIT(15) ,DATWFC(15) ,DATG(15) 
COMMON/C1/ME , M I , CHG , EPS 0 , K , HBAR , E6P , ESP , ACOEFF , PI , LS , 
*LP,SIGIA,SI GEA , DLTXCT , VD , JE , TI , TEM IT , AELEC , D I ST , G , NE , 
•SIGSP, SIGO i NNUC 
COMMON/C3 /KPS , KPC , KSP , KCP , NGS 
C0KM0N/C4/M0D , QP , SMI 
COMMON/C2/TE , DA , MUE 
COMMON / C5/REA , RE I ,R 
COMMON/ C6/SAHA 

OPEN ( UN I T= 7 , DEV I CE= ’ DSK ' , F I LE= ’ PL1 * ) 

C 

C PHYSICAL CONSTANTS 
C 

ME=5.686E-16 
MI=1 . 377E-10 
CHG=1 . 0 
EPSO=5.527E5 
K=8 .6176E-5 
HBAR=6. 5784E-16 
E6P=2 . 4625 
ESP=1.4319 
E0=E6P+ESP 
ACOEFF=3.57E7 
PI=3. 14159265 


> 


vj VT>K - '-V^. ,> , •, , 


LS=0.0 

LP=1.0 

SIGIA=1.2E-13 
SIGEA=4.E-14 
SIGSP=50.E-16 
SIG0=-1 • 01E-14 
C 

C SLOPE OF CROSS-SECTION FOR 6S-6P TRANSITION IN CESIUM (DLTXCT) 

C OBTAINED FROM NORCROSS & STONE: PP. 667, FIGURE 1, CURVE #4; REF. 29 
C NOLAN AND PHELPS, 1965. 

C 

DLTXCT=71E-16 

C 

C OTHER INPUTS 
C 

VD=1 . 0 
t:=i500. 

AELEC=1 . 

NNUC=1.E16 

C 

C READ AND STORE RUN PARAMATERS 
C 

ICOUNT=0 

100 ICOUNT=ICOUNT+l 

WRITE(5,475) ICOUNT 

475 FORMAT (2X, ’INPUT DATA:DIST,TEMIT,WFC,G' , 12X, 'CASE' ,21) 

READ (5, 476) DDIST (ICOUNT) ,DTEMIT( ICOUNT) ,DATWFC( ICOUNT) , 

$DATG( ICOUNT) 

476 FORMAT (4G) 

WRITE (5, 700) DDIST (ICOUNT) ,DTEMIT( ICOUNT) ,DATWFC( ICOUNT) , 

$ DATG ( I COUNT ) 

700 FORMAT (2X,4 (2X ,E12 . 5) ) 

IF (DTEMIT( ICOUNT) .NE.O.O)GO TO 100 
C 

C MAIN PROGRAM DO LOOP. 

C 

ICOUNT=ICOUNT-l 
DO 10 NX=1, ICOUNT 
DIST=DDIST (NX) 

TEM I T=DTEM I T ( NX ) 

WFC=DATWFC (NX) 

G=DATG (NX) 

JE=120.*TEMIT»»2«EXP(-WFC/(K*TEMIT) ) 

C 

C WRITE RELEVANT INPUT DATA AND HEADER TO OUTPUT FILE. 

C 

WRITE (7,480) NX 
WRITE(7,481)DIST 
WRITE (7 ,482 ) TEM IT 
WRITE(7,433)G 
WRITE (7, 484) WFC 
WRITE(7,485) JE 
WRITE(7,486)VD 
WRITE (7,487) 

480 FORMAT (1H1, ' INPUTS PL1 CASE',12,/) 

481 FORMAT (2X, 'DIST ='F5.3) 

482 FORMAT (2X, ’TEMIT='F8.2) 

483 FORMAT ( 2X , ' G -'F5.3) 

484 FORMAT ( 2X , ' WFC *'F5.2) 

485 FORMAT (2X, ' JE *'F6.2) 


H 


486 FORMAT (2X, ' VD ='F5.2) 

487 FORMAT ( ' • ,2X,2X, 'NE' ,8X, • XN1’ 

*,9X, ’TE',12X, 'SMI 

*8X, 'XJEN ' ,8X, ' RJ ' , 6X , ' R ' , 8X , ' REI ' ,8X, ’VK’ ,8X, 'VA' ,6X, 

* 'EFF' ,2X, 'PRAD/PIN* ,1X, ’SAHA') 

C 

C MAIN DO LOOP INDEXING NE 

C OPTION Is SCALING ORDERS OF MAGNITUDE 

C 

C DO 11 J=1 , 3 

C DO 12 1=1,9 

C NE= . 1E12* FLOAT ( I ) * 10 . * ‘FLOAT ( J ) 

C+ 

C OPTION 2: FIXED INCREMENT 

DO 11 1=1,30 
NE= . 05E15*FLOAT ( I ) 

C 

C »* "SOLVE THE AMBIPOLAR DIFFUSION EQUATION*** 

C 

c 

C FIRST FIND 'TMIN'. WHICH IS THE SOLUTION TO THE AMBIPOLAR DIFFUSION 
C EQUATION NEGLECTING RECOMBINATION. 

C THIS TEMPERATURE IS THE LOWER BOUND IN THE ROOTFINDING 
C FOR THE AMBIPOLAR DIFFUSION EQUATION WITH RECOMBINATION 
C (WHICH IS EVALUATED BY THE FUNCTION TDIFF) . 

C 

A=1500. 

B=5500. 

SAHA=0. 

CALL BISECT(TMIN,A,B, . 5,IFLAG) 

IF (IFLAG.LT. 2) GO TO 46 
WRITE (7 , 333) IFLAG 

333 FORMAT (2X,' ROOTFINDING FOR TMIN BOMBED’ ,2X, 21) 

GO TO 9 

46 CONTINUE 
TSAHA= (A+B) /2 . 

C 

C SOLVE THE AMBIPOLAR DIFFUSION EQUATION AND CALCULATE THE FOLLOWING 
C QUANTITIES IN THE PROCESS. 

C 

C 

C IF ROOTFINDING BOMBS THEN THE SOLUTION IS TOO CLOSE TO THE 
C LOWER LIMIT CALCULATED BY TMIN. CONDITIONS ARE CALCULATED 
C AT THE SAHA LIMIT. SAHA FLAG IS SET TO 1. 

C 

A=B 

B=5500. 

CALL BISECT (TDIFF ,A,B, .01, IFLAG) 

IF (IFLAG.LT. 2 )GO TO 47 
C WRITF (7 , 334) IFLAG 

C334 FOFJfAT(2X, 'ROOTFINDING FOR DIFF. EQN . BOMBED: 

C -DEFAULT TO SAHA VALUES’ ,2X, 21) 

3AHA=1. 

TE=TSAHA 
GO TO 212 

47 CONTINUE 
TE=(A+B)/2. 

ERROR-ABS (A-B) /2 . 


C 

C 


CALCULATION OF THE RATE CONSTANTS 


■.rjjsnsw- 


c 

212 CC=ME/(2.*K«TE) 

DD=ESP/(K*TE) 

VMSQ=2 . « ESP/ME 

C 

C FORMULATION FOR KSP BASED ON VARIABLE CROSS-SECTION. 

AA= (CC/PI )**(3./2.)*4.*PI*SIG0 

KSPA= (AA/ (2 . *CC) ) *EXP(-CC*VMSQ) • (VMSQ+1 ./CC) 

BB= (CC/PI ) • *1 . 5*2 . «PI • (ME*1 . E20) *DLTXCT 

KSPB=BB/ ( 2 . * CC ) • EXP ( -CC* VMSQ ) * ( ( VMSQ/1 . E2 0 ) * VMSQ+2 . * VMSQ 
*/(CC*l.E20)+2./(CC*l.E20*CC) ) 

KSP=KSPA+KSPB 

KPS=KSP»EXP (DD) * (2.*LS+1. )/(2.*LP+l. ) 

KCP=.46*EXP ( .5* (ALOG (8 . *TE/PI ) +ALOG (K) -ALOG (ME) ) -5 . • (ALOG (4 . *PI *TE 
• ) +ALOG (EPS0*K) ) ) 

KPC=KCP/ ( (2 . *LP+1 . ) • (2 . *PI • (HBAR/ME) « (HBAR/K) /TE) * * (3 ./2 . ) 
•*EXP(E6P/(K*TE))) 

R1=KSP/ ( KPC+KPS+G • ACOEFF/NE ) 

S=KSP-R1* (KPS+G-ACOEFF/NE) 

R0=NE«*2. * KCP/ ( KPS+KPC+G* ACOEFF/NE ) 

ALPHA=RO/NE« (NE*KPS+G*ACOEFF) /NE* *2 . 

NGS=NNUC-NE 
MFPI=1 ./ (NNUC*SIGIA) 

MFPE=1./(NNUC*SIGEA) 

MUI=2.*CHG*MFPI/(3.*(K*TI*PI«MI/2.)**.5) 

MUE=2 . *CHG*MFPE/ (3 . * (K*TE*PI*ME/2 . ) ** . 5) 

DA=MUI *K« (TE+TI)/CHG 
IF (SAHA. EQ. 0.0) GO TO 214 
QP=SQRT ( ALPHA/ (8.*DA))*NE*DIST 
IF(QP.LE.20.)GO TO 213 
C 

C IN THE CASE WHERE SMI (SM1=1-M0DULUS) IS TOO LOW TO BE CALCULATED BY 
C THE A.G.M. LADDER (QP>20) THEN SMI IS SET TO ZERO AND SAHA FLAG IS 2. 
C 

SAHA=2 . 

SM1=0. 

MOD=l . 

N6P~ (R0+R1*NGS) «AELEC*DIST 
GO TO 215 
C 

C ***CALCUALTE OUTPUT CURRENT BY THE ENERGY EQN . * * • 

C 

C 

C OUTPUT CURRENT CALCULATION (VARIABLE XJEN) : 

C N6P CALCULATED BY INTEGRATION OVER THE LENGTH. 

C THE INTEGRAND, N6P(X), IS CALCULATED BY THE EXTERNAL FUNCTION F. 

C 

213 SM1=16 . «EXP (-2 . *QP) 

M0D=1.-SM1 

214 CALL SIMPER (F, AREA, DIST) 

N6P=AREA 

215 PDIFF=4. *E0*AELEC*DA*QP*NE/DIST*1 . 602E-19 
PRAD=G*AC0EFF*ESP*N6P*1 .602E-19 
Q=PDIFF+PRAD 

C 

XJEN1=2 . * JE*K* (TE-TEMIT)/ (VD*CHG) 

XJEN2=Q/VD 

XJEN=XJEN1+XJEN2 

C 

C ***SOLVE THE TRANSPORT EQUATION* •• 


n o o 


c 

c 

C CALCULATE RJCRIT, THE UPPER LIMIT FOR ROOTFINDING, BASED ON THE CONDITION 
C THAT THE ELECTRON DENSITY AT THE ANODE XN1, MUST BE POSITIVE. 

C 

RJCRIT= (SIGIA/SIGEA) *SQRT ( (TI/TE) * (MI/ME) ) 

»* (CHG*AELEC)«DA«NE* (PI/DIST) *1.602E-19 
RJ=RJCRIT-RJCRIT* (.00001) 

I F (R J . GE . JE ) RJ=JE-JE* ( . 001 ) 

C 

C ROOTFINDING FOR OUTPUT CURRENT BASED ON THE TRANSPORT EQN. 

C 

C TRANSPORT EQUATION IS SOLVED BY EXTERNAL FUNCTION XJT. 

C THE VARABLE USED FOR THE OUTPUT CURRENT IS RJ. 

C 

A=RJ/100. 

B=RJ 

CALL BISECT (XJT, A, B, .001, IFLAG) 

I F ( I FLAG . LT . 2 ) GO TO 99 
WRITE (7, 717 )NE, IFLAG 

717 FORMAT ( ' ',X,E7.1,' ROOTFINDING XJT BOMBED' , 2X, • I FLAG= • ,21 ) 

GO TO 9 

99 CONTINUE 

RJ=(A+B)/2. 

ERROR=ABS(A-B)/2. 

CALCULATION OF SHEATH VOLTAGES 

VEBAR=SQRT(8.*K*TE/(PI*ME)) 

VIBAR=SQRT (8 . *K*TI/ (PI*MI) ) 

ALFA=SQRT (TE*2 . *PI/TI ) 

XNA= (4 . / (ALFA* VIBAR) ) «DA*2 . *QP*NE/DIST 

XNB= (SIGEA/SIGIA) *SQRT ( (TE/TI ) * (ME/MI ) ) «RJ/ (1 . 602E-19* AELEC) 

*• (4./ (ALFA* VIBAR) ) 

XN1=XNA-XNB 
XNO=XNA+XNB 

ARGVK=XN0*VEBAR*CHG*1 .602E-19/ (4. • (JE-RJ) ) 

VK=K « TE/CHG * ALOG ( ARGVK ) 

ARGVA=XN1* VEBAR*CHG* 1 . 602E-19/ (4 . *R J ) 

VA=K* TE/CHG* ALOG (ARGVA) 

C 

C WFC IS THE CATHODE WORK FUNCTION 
C WFA IS THE ANODE WORK FUNCTION. 

C 

WFA=WFC-VD 

C 

C CALCULATE SYSTEM EFFICIENCIES 
C 

EFF1=PRAD/ (PRAD+RJ* (WFA+2 . *K* TE/CHG) ) 

EFF2=PRAD/(RJ*VD) 

C 

C ESTABLISH ION CURRENT 
C 

JION=PDIFF/EO 

ratio=jion/rj 

c 

C WRITE calculated quantities 

c 

WRITE (7, 761 >NE,XN1,TE, SMI, XJEN,RJ,R,F:EI ,VK,VA,EFF1,EFF2, SAHA 
761 FORMAT ( • • ,X,E8. 3 , 2X,E9. 3, 3X, F10.4,4X. E9. 3 , 4X,F7 .2, 4X 


*,F7.2,3X,F5.3,3X,F7.4,3X,F7.3,3X,F7.3,3X,F5.3,3X,F6.3,3X,F2.0) 

9 CONTINUE 

12 CONTINUE 

11 CONTINUE 

10 CONTINUE 
VRITE(5,752) 

752 FORMAT ( 2X , ' OUTPUT WRITTEN TO FILE: PL1.DAT ') 

CLOSE (UNI T=5) 

STOP 

END 

FUNCTION TMIN(TE) 

C 

C CALCULATES IGNITED CONDITION BASED ON NO RECOMBINATION 
C 

REAL KSP , KPS , KSPA , KSPB , KCP , KPC , NGS , NE , MU I , MUE , MFP I , MFPE 
REAL MI,ME,HBAR,PI,LS,LP,K,NNUC,JE,JNET,MOD 
COMMON/C1/ME , M I , CHG , EPSO , K , HBAR , E6P , ESP , ACOEFF , PI , LS , 
*LP,SIGIA,SIGEA,DLTXCT, VD, JE,TI ,TEMIT,AELEC,DIST,G,NE, 

•SIGSP, SIGO.NNUC 
CC=ME/(2.*K*TE) 

DD=ESP/(K*TE) 

VMSQ=2.*ESP/ME 

AA= (CC/PI ) •* (3./2.)*4.*PI*SIG0 

KSPA= (AA/ (2 . *CC) ) *EXP (-CC* VMSQ) • (VMSQ+1 ./CC) 

BB= (CC/PI ) *»1 . 5*2. *PI* (ME*1 .E20) *DLTXCT 

KSPB=BB/ (2 . *CC) *EXP (-CC* VMSQ) « ( (VMSQ/1 .E20) * VMSQ+2 . *VMSQ 
*/(CC*l.E20)+2./(CC*l.E20*CC)) 

KSP=KSPA+KSPB 

KPS=KSP*EXP(DD)*(2.*LS+1.)/(2.*LP+1.) 

KCP=.46*EXP(.5* (ALOG(8.«TE/PI)+ALOG(K)-ALOG(ME) )-5.* (ALOG (4. *PI*TE 
* ) +ALOG (EPS0*K) ) ) 

KPC=KCP/ ( ( 2 . « LP+1 . ) • ( 2 . * P I • (HBAR/ME ) » (HBAR/K ) /TE ) * * ( 3 . /2 . ) 
••EXP(E6P/(K«TE))) 

R1=KSP/ ( KPC+KPS+G* ACOEFF/NE ) 

S=KSP-R1 • (KPS+G*ACOEFF/NE) 

RO=NE« * 2 . «KCP/ ( KPS+KPC+G* ACOEFF/NE ) 
ALPHA=R.7/NE»(NE*KPS+G*ACOEFF)/NE«»2. 

NGS=NNUC-NE 

TMIN=S*NGS-NE* ALPHA* NE 

RETURN 

END 

w 

FUNCTION TDIFF(TE) 

C 

C CALCULATES IGNITED CONDITION W/ RECOMBINATION ASSUMED. 

C 

REAL KSP , KPS , KSPA , KSPB , KCP , KPC , NGS , NE , M’JI , MUE , MFP I , MFPE 
REAL M I , ME , HBAR , PI , LS , LP , K , NNUC , JE , JNET , MOD 
COMMON/Cl/ME , M I , CHG , EPSO , K , HBAR , E6P , ESP , ACOEFF , PI , LS , 
«LP,SIGIA,SIGEA, DLTXCT, VD, JE,TJ r TEMIT, AELEC,DIST,G,NE, 
*SIGSP,SIGO,NNUC 
C0MM0N/C4/M0D , QP , SMI 
CC=ME/ (2 . *K*TE) 

DD=ESP/(K*TE) 

VMSQ=2.* ESP/ME 

AA= (CC/PI) •• (3./2.)*4.*PI*SIG0 

KSPA= (AA/ (2 . *CC) ) *EXP (-CC* VMSQ) • (VMSQ+1 . /CC) 

BB= (CC/PI) **1. 5*2. *PI * (ME*1.E20)*DLTXCT 

KSPB-BB/ (2 . «CC) «EXP (-CC* VMSQ) • ( (VMSQ/1 .E20) ‘VMSQ+2 . » VMSQ 


o o o 


*/(CC*l.E20)+2./(CC*l.E20*CC)) 

KSP=KSPA+KSPB 

KPS=KSP*EXP (DD) • (2 . «LS+1 . ) / (2 . *LP+1 . ) 

KCP= . 46*EXP ( . 5* (ALOG (8 . *TE/PI ) +ALOG (K) -ALOG (ME) ) -5 . • (ALOG (4 . «PI«TE 
• ) +ALOG (EPS0*K) ) ) 

KPC=KCP/ ( (2 . *LP+1 . ) • (2 . *PI • (HBAR/ME) • (HBAR/K) /TE) ** (3./2 . ) 

**EXP (E6P/ (K*TE) ) ) 

R1=KSP/ (KPC+KPS+G*ACOEFF/NE) 

S*KSP-R1* (KPS+G*ACOEFF/NE) 

R0=NE« *2 . *KCP/ ( KPS+KPC+G* ACOEFF/NE ) 

ALPHA=RO/NE* (NE*KPS+G*ACOEFF) /NE«*2 . 

NGS=NNUC-NE 

MFPI=1./(NNUC*SIGIA) 

MUI=2 . *CHG«MFPI/ (3.*(K*TI*PI«MI/2.)**.5) 

DA=MUI*K« (TE+TI)/CHG 
C 

C A.G.M. CALCULATION OF QUARTER PERIOD; QP 

C FOR MORE INFO SEE SECTION 16.4 IN THE HANDBOOK OF MATH. FUNCTIONS. 

C MODULUS (MOD) BASED ON RATIO OF NE (MAX) TO N (SAHA) . 

C 

MOD=l ./ (2. *S*NGS/ (NE*ALPHA*NE) -1. ) 

SM1=1 .-MOD 
A0=1. 

BO=SQRT(SMl) 

CO=SQRT (MOD) 

15 ANEW= (AO+BO) /2 . 

BNEW=SQRT (AO* BO ) 

CNEW= (AO-BO) /2 . 

ADIF=AO-ANEW 

IF(ADIF.LT. .00001) GO TO 16 

AO-ANEW 

BO=BNEW 

C0=CNEW 

GO TO 15 

16 CONTINUE 
QP=PI/(2.*ANEV) 

X1=S*NGS 

X2=4. *QP**2. *DA/DIST**2. 

X3= . 5*NE*ALPHA*NE 
TDI FF=X1-X2-X3 
26 RETURN 
END 
C 

FUNCTION XJT(RJ) 

CALCULATES THE TRANSPORT EQUATION 

REAL KSP , KPS , KCP , KPC , NGS , NE , MUI , MUE , MFPI , MFPE 
REAL MI,ME,HBAR,PI,LS,LP,K,NNUC, JE,JNET,MOD 
COMMON/C2/TE , DA , MUE 

COMMON/C1/ME , M I , CHG , EPSO , K , HBAR , E6P , ESP , ACOEFF , PI , LS , 
•LP,SIGIA,SIGEA,DLTXCT, VD, JE,TI ,TEMIT, AELEC,DIST,G,NE, 

•SIGSP, SIGO, NNUC 
COMMON/C4/MOD , QP , SMI 
COMMON/C5/REA , RE I , R 
COMMON/ C6/ SAHA 
VIBAR=SQRT(8.*K*TI/(PI*MI)) 

ALFA=SQRT (TE* 2 . *PI/TI ) 

XNA" (4./ (ALFA*VIBAR) ) *DA*2 . *QP*NE/DIST 

XNB= (SIGEA/SIGI A) *SQRT ( (TE/TI ) • (ME/MI ) ) *RJ/ (1 . 602E-19* AELEC) 


•• (4./ (ALFA*VIBAR) ) 

XN1=XNA-XNE 

XNO=XNA+XNB 

ARG1=NE*»2./ (XN0*XN1*4. ) 

TERM1-1 . /QP» ALOG ( ARG1 ) 

IF (SAHA. EQ. 0.0) GO TO 100 
TERM 2 =2 
GO TO 101 

100 TERM2=AL0G(16./SM1)/QP 

101 REA=DIST/ (2.»NE*MUE*1.602E-19) * (TERM1+TERM2) 

XLAM=12389 . «TE* *1 . 5/SQRT (NE/2 . ) 

REI=DIST«100/AELEC«ALOG (XLAM) / ( . 01 3085*TE**1 . 5) 
R=REA+REI 

C1=RJ*R 
C2=RJ/ (JE-RJ) 

C3=ALOG(C2) 

C4=C3*K*TE/CHG 
XJT= (C1+C4) -VD 
RETURN 
END 
C 

SUBROUT I NE B I SECT ( F , A , B , XTOL , I FLAG ) 

C 

C NO-FRILLS BISECTION ROOT-FINDING SUBROUTINE 
C 

IFLAG=0 

N*-l 

FA=F(A) 

FB=F (B) 

C VRITE(7,790)FA,FB 

C790 FORMAT ( 2X , 'ENDPOINTS: ' , 3X, ' FA= ' , E15 . 7, 7X, * FB= ' ,E15.7) 

IF (FA*F (B) .LE.O. ) GO TO 5 
IFLAG=2 
C 

C WRITE (7,601) A, B 

C601 FORMAT (43H F(X) IS OF SAME SIGN AT THE TWO ENDPOINTS 
C • 2F15.7) 

RETURN 

5 ERROR=ABS (B-A) 

6 ERROR=ERROR/2 . 

C 

C CHECK FOR SUFFICIENTLY SMALL INTERVAL 
C 

IF (ERROR. LE.XTOL) RETURN 
XM=(A+B)/2. 

C 

C CHECK FOR UNREASONABLE ERROR REQUIREMENTS 
C 

IF(XM+ERROR.EQ.XM) GO TO 20 
FM=F (XM) 

C 

C TEMPORARY PRINTOUT 
C 

N=N+1 

C 

C CHANGE TO NEW INTERVAL 
C 

IF(FA*FM.LE.O. ) GO TO 9 

A=XM 

FA=FM 


V: 


GO TO 6 
9 B=XM 
GO TO 6 
20 IFLAG=1 

KLTUKN 

END 

C 

FUNCTION F (X) 

C 

C CALCULATES N6p(X) = RO ♦ R1«NGS SO THAT IT CAN BE INTEGRATED BT 
C THE SUBROUTINE SIMPER. 

C 

COMMON /Cl /ME , M 1 , CHG , EPSO , K , HBAR , E6P , ESP , ACOEFF , PI , LS , 

•LP,SIGIA ; SIGEA,DLTXCT, VD, JE,TI ,TEMIT,AELF.C,DIST,G,NE, 

•SIGSP.SIGO , NNUC 
COMMON/C3/KPS , KPC ,KSP ,KCP , NGS 
COMMON/C4/MOD , QP , SMI 
DIMENSION A (20) ,B(20) ,C(20) ,PSI (20) 

REAL KPS , KPC , KSP , KCP , NGS , NE , MOD 
C 

C CALCULATE AND STORE A.G.M. SCALE 
C 

A(l)*l. 

B(1)=SQRT(SM1) 

C (1 ) =SQRT (MOD) 

N=2 

20 A(N)=(A(N-l)+B(N-l))/2. 

B(N) *=SQRT(A(N-1) »B(N-1) ) 

C(N)*(A(N-l)-B(N-I))/2. 

ADIF^A(N-1)-A(N) 

IF(ADIF.LT. .0001)GO TO 26 

N=N+1 

GO TO 20 

26 PSI (N)=2«* (N-l) »A (N) »2. »QP»X/DIST 

C 

C RECURSION CALCUIATION 
C 

DO 13 1*1, N 
NN=N-I+1 

TARG=C(NN)/A(NN) »SIN(PSI (NN) ) 

I F ( TARG . LT . 1 . 0 ) GO TO 29 
C WRITE (5,68) 

C WRITE(5,69)C(N),A(N) , PSI (NN) , TARG 

68 F0RMAT(2X, ’SUB:F(SAHA) BOMB: C(N) ,A(N) , PSI (NN) , TARG' ) 

69 FORMAT (2X,E15-9,4X,E15.9,4X,E15.9,4X,E15.9) 

TARG=1. 

29 PSI (NK-1) = (PSI (NN) +ASIN (TARG) ) /2 . 

13 CONTINUE 

SN=S1N (PSI (1) ) 

F* (SN» » 3 . «NE»KCP»NE* »2 .+NE , »SN»KSP»NGS) 

•/ (SN*NE* (KPS+KPC) +G*ACOEFF) 

RETURN 

END 

C 

SUBROUTINE SIMPER (FX, AREA, DIST) 

C 

C SIMPSONS RULE: NSEG MUST BE EVEN, FX MUST BE SYMMETRICAL ABOUT MIDPOINT. 
C 

REAL LSEG 
NSEG=20 


LSEG=DIST/FLOAT(NSEG) 

SUM=0 

X3=0. 

NSTEP=NSEG/4 
DO 39 1=1, NSTEP 
X1=X3 
FX1=FX(X3) 

X2=Xi+LSEG 

X3=X2+LSEG 

SUM=SUM+FXl+4. «FX (X2 ) +FX (X3) 
CONTINUE 

AREA=2 . *LSEG/3 . *SUM 

RETURN 

END 


